function F = solve_eqm_outer(param, const)
%% Initial Equilibrium

    % Start at the lowest quality
    q_indicator = sparse(param.Nf, param.Q);
    q_indicator(:,1) = 1;
    
    % Calculate the initial integrals
    ExportProbQ_init = ones(param.Nf,param.Q);
    integ = cal_integral(const, q_indicator, ExportProbQ_init);   
    
    
%% Solve the open economy model with quality choice

    % Setting
    tol = 1e-4;
    maxRepeat = 6;
    maxIter = 100;
    
    % Initialize
    eqm.Pi0 = ones(1,param.Q);
    eqm.Pi1 = ones(1,param.Q);  
    nq_new = cal_nq(const, q_indicator);
    nq_history = zeros(maxIter, param.Q);
%     figure;
    
    % Iteration
    err = 1;
    iter = 0;
    repeat = 0;  
    while err >= tol && iter < maxIter && repeat < maxRepeat      

        % update
        Pi0 = eqm.Pi0;
        Pi1 = eqm.Pi1;          
        nq = nq_new;
        iter = iter + 1;
                
        % solve equilibrium
        eqm = solve_eqm_inner(param, const, integ); 
        
        % update optimal quality 
        [q_opt_index, ~] = grid_search(param, const, eqm.Pi0, eqm.Pi1, eqm.ExportCutoffQ, eqm.ExportProbQ);        
        q_indicator = cal_indicator(param, q_opt_index); 

        % update the integral
        integ = cal_integral(const, q_indicator, eqm.ExportProbQ); 
        
        % update the endogenous industry structure
        nq_new = cal_nq(const, q_indicator);
        nq_history(iter,:) = nq_new;          
        
        % repeat count
        if iter > 2
            repeat = repeat + min(nq_new == nq_history(iter-2,:));    
        end                
        
        % check convergence 
        err1 = max(abs(nq_new-nq));
        err2 = max(abs(eqm.Pi0-Pi0));
        err3 = max(abs(eqm.Pi1-Pi1));        
        err = max([err1, err2, err3]);
        
%             % visualize progress
%             scatter(const.Qgrid, nq_new, 'Marker', '.')
%             xlabel('q')  
%             title(strcat('iter=',num2str(iter),', repeat=',num2str(repeat)))
%             drawnow limitrate
        
    end
    close
    

%% Store   
    
    eqm.iter_outer      = iter;
    eqm.err_outer       = err;
    eqm.repeat          = repeat;   
    
    eqm.q_index         = q_opt_index;
    eqm.q_indicator     = q_indicator;
    eqm.integ           = integ;  
    
    eqm.nq               = nq_new;
    eqm.first_nonzero    = find(eqm.nq, 1, 'first');
    eqm.last_nonzero     = find(eqm.nq, 1, 'last');  
    eqm.first_nonzero_q  = const.Qgrid(eqm.first_nonzero);
    eqm.last_nonzero_q   = const.Qgrid(eqm.last_nonzero);  
    
    eqm.density_all     = const.density_omega;  
    eqm.density_exp     = const.density_omega.*full(sum(eqm.ExportProbQ.*eqm.q_indicator,2));
    eqm.density_nonexp  = const.density_omega - eqm.density_exp;

    
%% Determine Quintile Group for each q

    % Quintileflag
    quintileflag = zeros(param.d, param.Q);
    quintilecutoff = zeros(1, param.d+1);
    nq_cum = cumsum(eqm.nq);
    for i = 1: param.d
        [~,quintilecutoff(i+1)] = min(abs(nq_cum - 0.2*i));
        quintileflag(i,quintilecutoff(i)+1:quintilecutoff(i+1)) = 1;            
    end
    quintileflag(1,(eqm.nq==0))=0;
    nd = sum(repmat(eqm.nq, param.d,1).*quintileflag,2)';  
        
    % Store
    eqm.nd                = nd;
    eqm.quintilecutoff    = quintilecutoff;
    eqm.quintileflag      = quintileflag;    
    
    
%% Indicators for Quintile Group for each firm type

    % Indicators for ex-ante quintile
    eqm.Qindicator = eqm.q_indicator*(eqm.quintileflag)'; %Nf by 5 matrix
    
    
%% Return

    F = eqm;

    
end